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A numerical model based on the finite-difference time-domain (FDTD) method is developed to 
simulate thermal noise in open cavities owing to output coupling. The absorbing boundary of the 
FDTD grid is treated as a blackbody, whose thermal radiation penetrates the cavity in the grid. 
The calculated amount of thermal noise in a one-dimensional dielectric cavity recovers the standard 
result of the quantum Langevin equation in the Markovian regime. Our FDTD simulation also 
demonstrates that in the non-Markovian regime the buildup of the intracavity noise field depends 
on the ratio of the cavity photon lifetime to the coherence time of thermal radiation. The advantage 
of our numerical method is that the thermal noise is introduced in the time domain without prior 
knowledge of cavity modes. 
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I. INTRODUCTION 

The finite-difference time-domain (FDTD) method |l| 
has been extensively used in solving Maxwell's equations 
for dynamic electromagnetic (EM) fields. The absorb- 
ing boundary condition based on the perfectly matched 
layer (PML) allows the simulation of open systems, 
e.g. leaky optical cavities, in any dimension. The 
incorporation of auxiliary differential equations, such 
as the rate equations for atomic populations [3| and 
the Maxwell-Bloch equations for the density-of-states of 
atoms 0, H, @] , has lead to comprehensive studies of light- 
matter interactions. Although the FDTD method has be- 
come a powerful tool in computational electrodynamics, 
it is applied mostly to classical or semiclassical problems. 
Recently, quantum fluctuations due to the spontaneous 
emission of atoms were introduced to the FDTD sim- 
ulation of microcavity lasers 0, HI- The light field in 
an open cavity also experiences quantum fluctuation be- 
cause of its coupling to external reservoirs. In this paper, 
we model the quantum noise for the cavity field as a clas- 
sical noise and incorporate it into the FDTD algorithm. 

There are two dissipation mechanisms for the cavity 
field: (i) intracavity absorption, (ii) output coupling. In 
the modal picture, widely used in quantum optical stud- 
ies, thermal noise is introduced so that the quantum oper- 
ator of a leaky cavity mode satisfies the commutation re- 
lation. Although thermal noise is quantitatively insignif- 
icant at optical frequencies, its proper treatment consti- 
tutes an essential part of the exact quantum-mechanical 
theory of lasers. Early laser theory introduces the ther- 
mal noise via a heatbath made up of loss oscillators or 
absorbing atoms [H, [l(|. It accounts for light absorp- 
tion inside the cavity. For a laser cavity whose loss only 
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comes from the output coupling, the thermal noise is 
attributed to the thermal radiation that penetrates the 
cavity through the coupling [lj, [ijj ■ Thus the amount 
of thermal noise depends on the mode decay rate, which 
must be known in order to solve the Langevin equation 
for the field operator. For open complex cavities, e.g. the 
ones made of random structures, the required informa- 
tion of modes is unknown a priori. Thus, it is desirable 
to be able to study the noise of a cavity field without 
prior knowledge of cavity modes. Additional problems 
with the modal picture are, (i) if the cavity is very leaky, 
the significant overlap of modes in frequency makes it dif- 
ficult to distinguish one mode from another; (ii) In the 
presence of nonlinearity, strictly speaking, the modes do 
not exist. In fact, one advantage of the FDTD method is 
the direct time-domain calculation of EM fields without 
prior knowledge of modes. The effective modal behavior 
is an emergent property that results from temporal eval- 
uation of the EM fields. We intend to introduce noise 
to the EM field in a way compatible with the FDTD 
method, namely, without invoking the modal picture. 
Our goal is to open a new approach for the study of 
quantum mechanical aspects of radiation in macroscopic 
systems with classical electrodynamics simulations. We 
believe our approach has the potential to permit rigorous 
theoretical investigations of noise in the area of quantum 
optics and of open systems such as chaotic open cavities. 
The dynamics of such systems are, in particular, very 
difficult to study using the standard frequency domain 
methods. 

In FDTD simulations, it is rather straightforward to 
introduce noise related to intracavity absorption. A fluc- 
tuating electric field can be added as a soft source at 
every grid point inside the cavity with its rms amplitude 
proportional to the local absorption coefficient [l3|, [l4[ • 
The output coupling, however, is not a local loss. The 
question is how to introduce thermal noise related to cav- 
ity leakage without knowing the leakage rate. In FDTD 
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simulations, light escaping from an open system is ab- 
sorbed by the absorbing boundary layer (ABL) which 
acts as the external reservoir. Since it absorbs all im- 
pinging fields, the ABL can be modeled as a blackbody. 
To remain in thermal equilibrium, the blackbody must 
radiate into the system. The blackbody radiation from 
the ABL propagates into the cavity and acts as noise to 
the cavity field. The amount of noise penetrating the cav- 
ity depends on the cavity openness or output coupling. 
We simulate the blackbody radiation from the ABL in 
the FDTD calculation. Our model is validated in the 
calculation of field noise in a one-dimensional dielectric 
cavity. In a good cavity whose lifetime r is much longer 
than the coherence time of thermal radiation r c , the av- 
erage amount of thermal noise in one cavity mode agrees 
to the solution of the quantum Langevin equation under 
the Markovian approximation. In addition to recovering 
the standard results, our simulations with various values 
of r and t c illustrate the transition from the Markovian 
regime to the non-Markovian regime, and demonstrate 
that the buildup of the intracavity noise field depends on 
the ratio of r c to r. This result is explained qualitatively 
by interference effect. 

The paper is organized as follows. Section |TT] outlines 
our numerical method. Possible numerical difficulties 
and problems are discussed. In Section IIIII we present 
the results of the FDTD calculations, including the black- 
body radiation in vacuum and noise penetration into a 
one-dimensional (ID) cavity. The transition from the 
Markovian regime to the non-Markovian regime is stud- 
ied. SectionHVlconsists of a discussion and interpretation 
of the results. An analytical expression is found which 
offers further insights to the amount of noise inside an 
open cavity. We end in Section [V] by summarizing all 
results and including a discussion of future applications 
for our method. 



II. NUMERICAL METHOD 

Our numerical model is based on the key insight that 
the ABL normally used to bound FDTD computational 
grids is in effect a blackbody which ideally absorbs all in- 
cident radiation. To stay in thermal equilibrium with 
temperature T, the blackbody must radiate into the 
system. To simulate the blackbody radiation, we sur- 
round the grid with a series of noise sources next to 
the grid/ABL interface. These soft sources radiate EM 
waves into the grid having spectral properties consis- 
tent with blackbody radiation. In this paper, we fo- 
cus on ID systems. The ID grid is discretized with a 
spatial step Ax and time step At. As shown in the in- 
set of Fig. [TJ two point sources are placed at the ex- 
tremities of the grid. Each source generates an electric 
field E s at every time step tj. Examples of the noise 
source of electric field E s (tj) are shown in Fig. [TJ A 
Fourier transform of the temporal correlation function of 
the electric field, (E s (ti)E s (t2)), gives the noise spec- 
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FIG. 1: (Color online) Noise source electric field E 3 (tj) gen- 
erated for T = 30, 000 K (black dots) and T = 50, 000 K 
(red crosses). The noise correlation time r c ~ 0.337 fs for 
T = 30, 000 K and r c w 0.203 fs for T = 50, 000 K. Ax = 1 
nm, M = 2 21 and r S i m = 7 ps. The inset is a schematic show- 
ing the noise sources placed next to the grid/ABL interface. 



trum D(u),T). If E s (tj) is uncorrelated in time, i.e., 
(E s (ti)E s (t 2 )) oc S(t 2 - ti), D(u,T) is the white-noise 
spectrum. This is incorrect as D{u),T) should be equal 
to the energy density of the blackbody radiation , 
which in one dimension is 

^ ' ^ 7rc (exp(/iu; / kT) — l) ^ ^ 

D{lo,T) for two different temperatures is plotted in 
Fig. [21 For computational convenience, we extend the 
range of lj from (0, oo) to (—00,00). Since the elec- 
tric field in the FDTD simulation is a real number, 
D(—lo, T) must be equal to D(cu, T) for u> > 0. Therefore, 
D(uj,T) = D(\uj\,T). We normalize D{uj,T) as 

so that f x oo D n (\uj\,T)duj = 2ir. 

The temporal correlation function for the source elec- 
tric field is given by 

r2 poo 

(E s (h)E s (t 2 )) = — / dwA.CM.TJeMfe-ti), (3) 

27T J -co 

where S is the rms amplitude of the noise field whose 
value is to be determined later. For the thermal noise, 
the field correlation function is given specifically by 

IS 2 

(E.itjE.fo)) = — [C(2, 1 - i(t a - h)kT/h) 

+£(2,l + t(i 2 -ti)fcT/fi)], (4) 

where the ^-function is given as 

00 

C( S ,a)=5> + a)- s . (5) 

k=0 
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FIG. 2: (Color online) FDTD-calculated energy density of 
blackbody radiation propagating in ID vacuum versus fre- 
quency u) for temperatures T = 30, 000 K (lower) and T = 
50,000 K (upper). The inset shows the energy density for 
temperature T = 30, 000 K at higher frequencies. The data 
are obtained by averaging over 2000 calculations with the res- 
olutions Ax = 10 nm (red crosses) and Ax = 1 nm (black 
dots). The source spectra D(u>,T) are also plotted as solid 
lines on top of the numerical spectra. 
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FIG. 3: (Color online) Temporal correlation function, 
(E s {ti)E s (t2)) vs. <2 — ti, for the noise electric field at 
T = 30, 000 K (black circles) and T — 50, 000 K (red crosses). 
The noise correlation times are r c w 0.337 fs for T = 30, 000 
K and r c « 0.203 fs for T = 50, 000 K. Ax = 1 nm, M = 2 21 , 
and T S im = 7 ps. The lines represent (E s (ti)E s (t2)} given by 
the analytical expression in Eq. [4] for T — 30, 000 K (black) 
and T = 50, 000 K (red). Every 5th data point is taken from 
the numerical data in order to better show the agreement with 
the analytical solution. 



The temporal correlation function of thermal radiation 
is plotted in Fig. [3] 

We employ a quick and straightforward way of generat- 
ing random numbers for E s (tj) so that Eq. [3] is satisfied. 
Freilikher et al. have developed such a method in the 
context of creating random surfaces with specific height 



correlations [l7|. The end result takes advantage of the 
fast Fourier transform (FFT) which we use to generate 
the source electric field: 
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where 2M is the total number of time steps, T S i m = 
2M At is the total simulation time and uji = 2itl/T S i m . 
Mi and Ni are independent Gaussian random numbers 
with zero mean and a variance of one. Their symmetry 
properties are Mi — M-i and Ni = — iV_j. These Gaus- 
sian random numbers can be generated by the Marsaglia 
and Bray modification of the Box-Miiller Transformation 
[HI], a very fast and reliable method assuming the uni- 
formly distributed random number generator is quick and 
robust. 

The electric field sources generate both electric and 
magnetic fields, which propagate into the grid. E(x, lu) 
and H(x,u>) are obtained by the discrete Fourier trans- 
form (DFT) of E(x,t) and H(x,t). Since both E(x,t) 
and H(x,t) are real numbers, E(x,u>) — E(x 1 —uj) and 
H(x,u>) = H(x,—u>). The EM energy density at fre- 
quency u> shall include E(x,u>), E(x,—oS), H{x,lu) and 
H(x,—u>). If the grid is vacuum, the steady-state en- 
ergy density at every position x should be equal to the 
blackbody radiation density. The rms amplitude 5 of the 
source field E s is determined by 
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When setting the parameters in the FDTD simula- 
tion, we must taken into consideration the characteris- 
tics of thermal noise. The temporal correlation time or 
coherence time t c of thermal noise is defined as the full 
width at half maximum (FWHM) of the temporal field 
correlation function. If the time step At is close to r c , 
E s exhibits a sudden jump at each time step. The ID 
FDTD algorithm cannot accurately propagate such step- 
like pulses (with sharp rising edge) if the Courant factor 
S = cAt/ Ax is set at a typical value S < 1. The pulse 
shape is distorted with fringes corresponding to both re- 
tarded propagation and superluminal response This 
occurs because the higher frequencies from the step dis- 
continuity are being inadequately sampled and because of 
numerical dispersion arising from the method of obtain- 
ing the spatial derivatives for E and H . To avoid such 
problems, we use S = 1 which eliminates the numeri- 
cal dispersion artifact [l[. Furthermore, we set At <C r c 
which provides a dense temporal sampling relative to the 
correlation/coherence time of the thermal noise. 

To obtain an accurate noise spectrum with the DFT, 
both the frequency and temporal resolutions must be cho- 
sen carefully. The two problems affecting the reliability 
of the DFT are aliasing and leakage due to the use of a fi- 
nite simulation time [191 ] . The solution to these problems 
is to increase the number of time steps 2M and decrease 
the time step value At. This takes the DFT closer to 



4 



a perfect analytical Fourier transform, but run-time and 
memory limitations must be considered as well. Tak- 
ing advantage of the FFT algorithm significantly reduces 
both noise generation time and spectral analysis time. 

Although the thermal noise spectrum can be very 
broad, only noise within a certain frequency range is rel- 
evant to a specific problem. Let ui m in and u> max denote 
the lower and upper limits of the frequency range of in- 
terest, and Auj the frequency resolution needed within 
this range. To guarantee the accuracy of the noise simu- 
lation in uj min < lu < Lo max , the total running time T sim 
must exceed 2u juj rn i n and 2u j Au>. The time step At has 
an additional requirement, At < ir/u> max . 



III. SIMULATION RESULTS 

A. Blackbody radiation in vacuum 

We first test the noise sources in a ID FDTD system 
composed entirely of vacuum. Two sets of independent 
noise signals E s (tj) are generated via Eq. O One set is 
added as a soft source one grid cell away from the left 
absorbing boundary; the other one cell from the right 
absorbing boundary. Both have equal rms amplitude S 
so that the average EM flux to the left equals that to the 
right at any position x in the grid. Since the system is 
one dimensional, the EM flux at any distance away from 
the source has the same magnitude. The value of 6 shall 
be adjusted so that Eq.[7Jis satisfied. For the EM energy 
density radiated by one source to equal D(\w\,T), we set 
S to 

After the noise fields in the grid reach steady state, the 
noise spectrum at any grid point is obtained by a DFT. 
We verify that the spectrum of EM energy density at 
any point in the grid is identical to that at the source. 
This means there is no distortion of the noise spectrum 
by the propagation of noise fields in vacuum. The two 
point sources at the grid/ABL interface radiate into both 
the grid and ABL. Since the two sources are uncorrelated 
with each other, their energy densities, instead of their 
field amplitudes, add in the grid. Thus no further modi- 
fication of 5 from that given by Eq. [5] is needed to satisfy 
Eq. \7\ It is numerically confirmed that 8 does not depend 
on the total length of the system. 

Examples of the noise source of electric field E s (tj ) at 
T = 30, 000 K and T = 50, 000 K are shown in Fig. Q] 
Air = 1 nm, and At -C r c . The frequency range of in- 
terest is set as uj m i n — 2 x 10 15 Hz, uj max — 2.5 x 10 16 
Hz, and the frequency resolution Auj = 1 x 10 12 Hz. 
^From the condition At < ir/u> max , Ax shall be less than 
37 nm. Figure O compares the FDTD calculated energy 
density to that of thermal radiation density D(co,T). Us- 
ing Ax = 10 nm creates a slight discrepancy at high 



frequencies; e. g. at u> < 1 x 10 16 Hz the mean er- 
ror > 2.5%. To reduce the error to below 2.5% at 
u m ax = 2.5 x 10 16 Hz, we refine the resolution. Using 
Ax = 4 nm changes the error at uj max to 1.6%. If the to- 
tal time step 2M = 2 21 is fixed, the decrease of At leads 
to a reduction of T S i m = 2MAt, which increases 27r/r S j m 
to 2 x 10 11 Hz. We must check that 27r/r S i m < ui m i n 
and 27r/r S i m < Auj are still satisfied. With Ax = 1 nm, 
the error at 2.5 x 10 16 Hz is further reduced to < 0.1%. 
27r/r s i TO increases to 9 x 10 11 Hz, which is still below the 
set values of u) m i n and Auj. Therefore, using the value of 
S in Eq. [5] and carefully choosing the numerical resolu- 
tions yield the blackbody spectrum at every point in the 
grid within the frequency range of interest. 

Figure [3] compares the FDTD calculated temporal cor- 
relation function of the electric field to that of Eq. 0] 
at T = 30,000 K and 50,000 K. With increasing tem- 
perature, the coherence time t c reduces quickly. The 
quantitative dependence of r c on T is found to be r c ss 
1.32fi,/ kT. This l/T dependence does not change for a di- 
mensionality higher than one; only the prefactor changes 
(20| . As the correlation time r c decreases, the time step 
At shall be reduced to maintain the temporal resolution 
of the correlation function. The subsequent reduction 
of total running time T s i m — 2MAt does not affect the 
numerical accuracy, as long as the total number of time 
steps 2M is fixed. A decrease of 2M would result in 
an increased mean-square error in the correlation func- 
tion due to less sampling. As shown in Fig. [3J the good 
agreement of the FDTD calculated temporal correlation 
function to that of blackbody radiation given by Eq. [T] 
confirms that introducing noise sources with the charac- 
teristics of blackbody radiation at the FDTD absorbing 
boundary effectively simulates thermal noise in vacuum. 



B. Thermal noise in a ID cavity 

Next we calculate the thermal noise in a dielectric slab 
of length L and refractive index n > 1 . This slab consti- 
tutes an open cavity in that electromagnetic field leakage 
occurs from both surfaces of the slab into an exterior re- 
gion. A schematic of the ID open cavity is shown in 
the inset of Fig. 4(b). The cavity mode frequency is 
u) m = mirc/nL, where m is an integer and c is the speed 
of light in vacuum. The frequency spacing of adjacent 
modes is du> = uj m — uj m -i = irc/nL, which is indepen- 
dent of m. Ignoring intracavity absorption, the decay 
of cavity photons is caused only by their escape from 
the cavity. All the cavity modes have roughly the same 
decay time r = 1/fejC, where ki = — In (r 2 ) /2nL, and 
r = (1 — n)/(l + n) is the reflection coefficient at the 
boundary of the dielectric slab. The mode linewidth is 
Suj = 2/t. We simulate only good cavities whose modes 
are well separated in frequency, namely, 5uj < duj. Since 
Suj oc 1/L, the ratio Suj/duj is independent of L, and is 
only a function of n. 

The Langevin equation for the annihilation operator 
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o-m(t) of photons in the m-th cavity mode is 

77 — = a m (t) + F m {t), 

at t 

where F m (t) is the Langevin force. If the noise cor- 
relation time r c <C t, F m (t) can be considered 5- 
corrclatcd in time. The Markovian approximation gives 

(Fh(t)F m (tf)) = D F S(t - t'). According to the fluc- 
tuation dissipation theorem, Dp — (l/r)nr(a; m ), where 
nT{io m ) = (e huJm l kT — is the number of thermal pho- 
tons in a vacuum mode of frequency co m at temperature 
T 0. 

From Eq. [5J the average photon number in one cavity 
mode (n m (t)) = (al n (t)a m (t)') satisfies 

d 2 2 

37 (n m (t)) = -- (n m {t)) + ~n T {u m ). (10) 
at t t 

At steady state, (n m ) = nr(w m ) in each cavity mode. 
The number of thermal photons is determined by the 
Bose-Einstein distribution nr(w m ). (h m ) is independent 
of the cavity mode decay rate because the amount of 
thermal fluctuation entering the cavity increases by the 
same amount as the intracavity energy decay rate. 

In our FDTD simulations, we verify that when r 3> t c 
the number of thermal photons in one cavity mode is 
equal to riT{uj m )- Since there is neither a driving field 
(e.g. a pumping field) nor excited atoms in the cavity, 
the EM energy stored in one cavity mode comes entirely 
from the thermal radiation of the ABL which is coupled 
into that particular mode. The steady-state number of 
photons in the m-th cavity mode is obtained from the 
FDTD calculation of intracavity EM energy within the 
frequency range w m _ 1/2 < w < lu„ 1+1/2 , where w m±1 / 2 = 
(m±l/2)irc/nL. 

I f«m-H/2 

n m = (h m ) = / du 



-1/2 



dx(^e\E(x,u)\ 2 + ^ \H(x,u)\ 2 ) (11) 

In our simulation, the temperature of the thermal sources 
at the ABL is T = 30,000 K. The coherence time of 
thermal radiation is r c = 0.337 fs. The refractive index 
of the dielectric slab is n — 6, and the length is L = 
2400 nm. The cavity lifetime t = 143 fs, is much longer 
than t c . The reason to choose a relatively large value 
of n is to have Su> < du> so that the cavity modes are 
separated in frequency. Care must be taken in setting 
the grid resolution because the intracavity wavelength is 
reduced to X/n. To maintain the spatial resolution, Aa; is 
decreased to meet Ax <C X/n. In our FDTD simulation, 
Ax = 1 nm and 2M = 2 21 . After the intracavity EM 
field reaches the steady state, we calculate the average 
thermal energy density inside the cavity 



Figure BJa) shows the intracavity noise spectrum U(w), 
which features peaks at the cavity resonant frequencies 
(9) oj m - Because n > 1, EM energy is also stored in the di- 
electric slab at frequencies away from cavity resonances. 
For example, U{lo = uJ m ±\/2) is higher than that in vac- 
uum by a factor of 2n 2 /(n 2 + 1). Thus the entire noise 
spectrum lies on top of the vacuum blackbody radiation 
spectrum, as confirmed in Fig. f?] The number of ther- 
mal photons in a cavity mode is calculated via Eq. [Tl] 
and plotted in Fig. [5] The modal photon number n m is 
equal to riT^LOm) with a mean error less than 0.1%. This 
result agrees with the steady-state solution of Eq. [TTJ It 
confirms that our numerical model of thermal noise in an 
open cavity is consistent with the prediction of quantum 
mechanical theory. Note that the modal photon numbers 
in Fig. 5 are time averaged values. Their values being 
much less than unity can be interpreted in a quantum 
mechanical picture as that most of the time there is no 
photon in the cavity mode. 



The above calculation is done in the Markovian regime. 
Next we move to the non-Markovian regime by reduc- 
ing t. The refractive index is kept at n = 6, while 
the cavity length L is reduced. This is a simple way 
of increasing the mode linewidth 5u> while keeping the 
modes separated in frequency, i.e. keeping Sui/duj con- 
stant. If r is decreased to less than r c , At shall be re- 
duced to keep At <C r. Meanwhile, the increase of the 
mode linewidth and mode spacing allows low frequency 
resolution, namely, an increase of Aui. For example, at 
L = 20 nm, we set Ax = 0.1 nm, Auj = 9 x 10 12 Hz 
and 2Af = 2 21 . Figure 0Jb) shows the intracavity noise 
spectrum U(lu) in this regime, which also features peaks 
at the cavity resonant frequency u> m . Figure [5] shows 
the FDTD-calculated value of n m as L decreases grad- 
ually from 2400 nm to 20 nm. When r approaches r c , 
n m is no longer independent of t, but starts increasing 
from nr(w m ). This means the number of thermal pho- 
tons that are captured by a cavity mode increases with 
the decrease of r. As the coherence time of thermal ra- 
diation impinging onto the cavity approaches the cavity 
photon lifetime, the constructive interference of the ther- 
mal field is improved inside the cavity, leading to a larger 
buildup of intracavity energy. 



We also investigate a different situation where r is 
fixed and r c is varied. By reducing the temperature T, 
the coherence time of thermal radiation r c is increased. 
Meanwhile, the energy density of thermal radiation is 
decreased. As shown in Fig. [51 the number of thermal 
photons n m in a cavity mode is reduced. This can be 
easily understood as there are fewer thermal photons in- 
cident onto the cavity at lower T. Nevertheless, n m is 
larger than n-r(w m ) at the same T. This is because the 
longer coherence time of the thermal field results in bet- 
(12) ter constructive interference inside the cavity. 
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FIG. 5: (Color online) The number of thermal photons in in- 
dividual cavity modes n m , calculated via the FDTD method, 
for a slab cavity with n — 6. The cavity length L is varied to 
change r. The impinging blackbody radiation has T = 30, 000 
K and r c = 0.337 fs. The values of r c /r are 0.0024, 0.29, 0.43, 
and 0.56. Lines are drawn to connect the data points at the 
mode frequencies uj m — ncm/nL to illustrate its frequency 
dependence. For t c <t (only every 5th mode for r = 143 fs 
is shown to improve the visibility), the photon number n m co- 
incides with the Bose-Einstein distribution nr. When r c ~ r, 
n m deviates from nr. 



FIG. 4: (Color online) Spatially-averaged EM energy density 
U(u>), calculated by the FDTD method, versus frequency uj in 
a dielectric slab cavity with refractive index n — 6 and length 
L = 2400 nm (a) and L = 20 nm (b). The vertical black 
dashed lines mark the frequencies of cavity modes u m . The 
spectrum of impinging blackbody radiation D(u),T) is also 
plotted (red solid line). In (a) the cavity decay time r = 143 
ps, much longer than r c . In (b), r = 1.19 fs, comparable to 



IV. DISCUSSION 

To gain a better understanding of our FDTD 
simulation results in the non-Markovian regime, we 
analytically examine the effect of noise correlation 
time t c on the amount of thermal noise inside an 
open cavity. The ratio of the intracavity EM 
energy at frequency lu to the energy density of 
the thermal source outside the cavity is W{u>) = 

{j L dx [±e\E(x,uj)\ 2 + i/iol^MI 2 ]} /D(u,T). For a 

dielectric slab of refractive index n and length L, we ob- 
tain the exp ression for W(u>) using the transfer matrix 
method EflL 



W(u) = 



2ujnL{l + n 2 ) / c + (n 2 - l)sin(2wni/c) 
1 + 6n 2 + n 4 — (n 2 — l) 2 cos(2umL / c) 
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FIG. 6: (Color online) The number of thermal photons in in- 
dividual cavity modes n m , calculated via the FDTD method, 
for a dielectric slab cavity with n — 6 and L — 20 nm. The 
cavity decay time r = 1.19 fs. The temperature of blackbody 
radiation is varied to change r c . The values of r c /r are 0.29 
(T = 30, 000 K), 0.43 (T = 20, 000 K), and 0.56 (T = 15, 000 
K). Black dashed lines are drawn to connect the data points 
at the mode frequencies U) m ~ ircm/nL to illustrate its fre- 
quency dependence. For comparison, the Bose-Einstein dis- 
tribution tit(w) is also plotted (red solid lines) for the same 
temperatures. 



(13) 
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It can be used to calculate the ratio B m {r c ,T) 
n m /n T (uj m ), as 



a + l/2 



(1ujW(uj)D(lu,T) 



/huj Tl 



(14) 



In the Markovian regime t c <t and D(u>,T) is nearly 
constant over the frequency interval of one cavity mode 
so 



S m (r c ,T) = 



D(gj m ,T) 
hw m n T (u) m ) 



w m + l/2 



dio W{w) 



1 

nc 



1 + 1/2 



dwW{oj) 



(15) 



-1/2 



We input the same parameters as the FDTD simulation: 
n = 6, L = 2400 nm, and r = 143 fs. As r c approaches 
zero, the integration of W(u>) from w m -i/2 to w m +i/2 
gives a value close to ire. Thus, as shown in the inset 
of Fig. [3 B m (r c , t) ps 1 for r c /r < 1. The deviation of 
B m {r c , t) from one is greater for smaller m. One possible 
reason is that the condition 5u> <C u m no longer holds for 
small m and there is a large uncertainty in defining the 
frequency of a cavity mode whose linewidth is compara- 
ble to its center frequency. In other words, the calculation 
of n m using Eq. 1141 becomes questionable. 

In the non-Mar kovian regime r c > r, D(oj,T) is not 
constant over the frequency range of a cavity mode any- 
more, thus it cannot be taken out of the integral in Eq.[T4l 
The behavior of B m (T c , r) in this regime is shown in the 
main panel of Fig. [7] As r c approaches r, B m (r c , t) no 
longer stays near one but increases with t c /t. This result 
is consistent with that of the FDTD simulation presented 
in the previous section. On one hand, if r is fixed and 
t c is increased by decreasing the temperature T, the ab- 
solute number of thermal photons in a cavity mode n m 
decreases, but its ratio to the number of thermal pho- 
tons in a vacuum mode n-r(wm) increases. On the other 
hand, if r c is fixed and r is decreased by shortening the 
cavity length L, both n m and n m /nT{w m ) increase. The 
departure of n m from riT(uj m ) is a direct consequence of 
the breakdown of the Markovian approximation. When 
the coherence time of thermal radiation is comparable to 
the cavity decay time, the Langevin force F m {t) in Eq. [S] 
is no longer (5-correlated in time, and Eq. [10] is invalid. 



V. CONCLUSION 

We have calculated the fluctuations of EM fields in 
open cavities due to output coupling with the FDTD 
method. The fluctuation dissipation theorem dictates 
that the cavity field dissipation by leakage be accompa- 
nied by thermal noise, which is simulated here by classical 
electrodynamics. The absorbing boundary of the FDTD 
grid is treated as a blackbody, which radiates into the 
grid. We have synthesized the noise sources whose spec- 
trum is equal to that of blackbody radiation. Careful 
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FIG. 7: (Color online) B m (r c ,r) = n m /nT(^m) vs. the ratio 
of noise correlation time to cavity decay time r c /r. B m (r c , r) 
depends only on the ratio r c /r and not on r c nor r individ- 
ually. This plot was generated by fixing r and varying r c . 
Cavity modes of higher m have a larger value of B m (r c ,r) 
whether in the regime t c ~ r (main panel) or r c r (inset). 



selection of numerical parameters in the FDTD simula- 
tion avoids the distortion of the noise spectrum by wave 
propagation in the ID grid. It is numerically confirmed 
that the noise fields propagating in vacuum retain the 
blackbody spectrum and temporal correlation function. 
When a cavity is placed in the grid, the thermal radi- 
ation is coupled into the cavity and contributes to the 
thermal noise for the cavity field. We calculate the ther- 
mal noise in a ID dielectric slab cavity. In the Marko- 
vian regime where the cavity photon lifetime r is much 
longer than the coherence time of thermal radiation t c , 
the FDTD-calculated amount of thermal noise in a cavity 
mode agrees with that given by the quantum Langevin 
equation. This validates our numerical model of ther- 
mal noise which originates from cavity openness or out- 
put coupling. Our FDTD simulation also demonstrates 
that in the non-Markovian regime the steady-state num- 
ber of thermal photons in a cavity mode exceeds that in 
a vacuum mode. This is attributed to the constructive 
interference of the thermal field inside the cavity. 

The advantage of our numerical model is that the ther- 
mal noise is added in the time domain without the knowl- 
edge of cavity modes. It can be applied to simulate com- 
plex open systems whose modes are not known prior to 
the FDTD calculations. Our approach is especially use- 
ful for very leaky cavities whose modes overlap strongly 
in frequency, as the thermal noise related to the cav- 
ity leakage is introduced naturally without distinguishing 
the modes. Therefore, we believe the method developed 
here can be applied to a whole range of quantum optics 
problems. Although in this paper the FDTD calculation 
of thermal noise is performed on ID systems, the exten- 
sion to 2D and 3D systems is straightforward. We com- 
ment that our approach does not apply to the simulation 
of zero-point fluctuation which has a different physical 
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origin from the thermal noise. However, our numerical 
method can be used to study the dynamics of EM fields 
which are excited by arbitrarily correlated noise sources. 
One potential application is noise radar [12, [Hj|. The 



propagation, reflection and scattering of ultra-wideband 
signals utilized by noise radar can be easily simulated 
using the method developed here. 
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